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In thermoacoustic tomography an object is irradiated by a short electromag- 
netic pulse and the absorbed energy causes a thermoelastic expansion. This ex- 
pansion leads to a pressure wave propagating through the object. The goal of 
thermoacoustic tomography is the recovery of the initial pressure inside the object 
from measurements of the pressure wave made on a surface surrounding the ob- 
ject. The time reversal method can be used for approximating the initial pressure 
when the sound speed inside the object is variable (non-trapping as well as trap- 
ping). This article presents error estimates for the time reversal method in the cases 
of variable, non-trapping sound speeds. Numerical examples for non-trapping as 
well as for trapping sound speeds are provided. 

1 Introduction 

Thermoacoustic tomography is a hybrid medical imaging technique characterized by 
high resolution and contrast. A short electromagnetic pulse is used to irradiate the bi- 
ological object of interest. A part of the electromagnetic energy is absorbed by the 
tissue, which causes a thermoelastic expansion. This leads to a pressure wave propa- 
gating through the object, which is measured by ultrasound transducers located on an 
observation surface S, usually surrounding the object. The collected information is 
used to reconstruct the initial pressure, which is roughly proportional to the absorbed 
energy [1,2]. The good contrast in the resulting images is due to electromagnetic 
energy being preferentially absorbed by cancerous cells, while ultrasound provides a 
submilimeter resolution [2,3]. 

Let 5 be a smooth observation surface surrounding the object and let B be the 
domain bounded by S. In applications only space dimensions n ^ 2,3 are of interest, 
but the analysis developed in this paper is carried out for an arbitrary dimension n. We 
assume that the speed c{x) of ultrasound in the tissue is smooth and strictly positive 
c{x) > c > 0, and that c{x) = 1 for large values of |a;|. Then the pressure p{x, t) at 
location x and time t satisfies the wave equation (e.g. [1,4-6]): 
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Here pa denotes the second time derivative of p{x, t), g{y, t) is the measured data, i.e. 
the value of the pressure at time t measured at transducer's location y £ S, and f{x) is 
the initial pressure, which is to be reconstructed. The function f{x) is assumed to have 
a compact support in B. One thus faces the problem of inverting the mapping R: f ^ 
g from the initial pressure to the measured data. For a more detailed description of the 
mathematical methods used in thermoacoustic tomography, one can refer to [6-1 1] and 
the references therein. There are various types of reconstruction procedures for closed 
observation surfaces, e.g. filtered backprojection formulas, eigenfunction expansion 
methods and time reversal method. A comparison of the three and a discussion of their 
advantages and limitations can be found in [12]. It was argued there that time reversal 
is the most versatile and easy to implement among these. Here we consider the errors 
involved in the time reversal method and provide estimates that justify its validity. 
It should be noted that, although we are interested in applications to thermoacoustic 
tomography, the time reversal method for the wave equation has been used in other 
applications [13-15]. Our results apply in all these cases. 

In Section 2 a non-trapping condition on the sound speed is recalled and the time 
reversal method is described. In Section 3 the main result of this artice is formulated 
and proved. Numerical examples are provided in Section 4. Section 5 contains some 
concluding remarks. This is followed by the acknowledgment section. 

Additional discussions and implementations of the time reversal method can be 
foundin[12,16,17]. 

2 Non-trapping condition and time reversal method 

In this section we give a more detailed description of the time reversal method and 
recall the notion of a non-trapping sound speed. 

2.1 Non-trapping condition 

We will be interested in the initial value problem: 



where c{x) > c > 0, c{x) € C°°(R"), and c{x) - 1, as well as fi{x) and f2{x), 
has compact support. Consider the following Hamiltonian system in 2n real variables 

(a;, C) with Hamiltonian H = ^ \^\^: 



The solutions of this system are called bichamcteristics and their projections into the 
a;-space are called rays. 

Deiuiition 1. We say that the non-trapping condition holds, if all rays (with Co 7^ 
tend to infinity when t ^ oo . 



Ptt = c^{x)Ap, t>0, xeW^ 
p{x,0) = fi{x), ptix,0) = f2{x), 



(2) 
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It is known that the singularities of the solution of (2) propagate along bicharacter- 
istics (e.g. [18-20]). Therefore, due to the non-trapping condition imposed on c{x), 
it follows that for any distributions /i , /2 with compact supports, the singularities of 
the solution move away to infinity as t ^ oo. Then, for a sufficiently large t, the so- 
lution p{x, t) is infinitely differentiable inside the unit ball B. Moreover, the solution 
decreases inside B and the following estimates hold: 

Theorem 2. ( e.g. [20, 21 ]) Under the conditions imposed on the sound speed c{x) and 
f{x) and for any bounded domain B, the solution of (2) satisfies the estimates 

f)k+\m\ 

^ < Crjkit) (||/i||l2 + II/2IU2) ,x€B,t>To, (4) 



for any multi-index (k,m), where rjk{t) = t~^^^~'^ for even n, and rjk{t) = e~^* for 
oddn. Here S is a positive constant depending only on c{x), Tq depends on the domain 
B, and the C depends on B and the multi-index (k, m). 

We will make use of the above theorem in estimating the error of reconstruction of 
f{x) when using time reversal. Note that in the case of trapping sound speed the local 
energy of the solution of (2) still decreases in any compact domain, but, in general, 
there is no uniform local energy decay estimate [22, 23]. 



2.2 The time reversal method 

The Huy gens' principle states that in odd dimensions and when the sound speed is 
constant, for any initial source with bounded support and for any bounded domain, 
there is a moment in time when the wave leaves the domain (e.g., [24]). Thus, given an 
initial pressure f{x) with a bounded support, there is a time T when the wave inside 
the domain B, bounded by the observation surface S, vanishes for all t > T. Then, 
to reconstruct f{x) we can simply "rewind" the solution, i.e. we can solve the wave 
equation backwards in time inside B with zero initial conditions att = T and boundary 
conditions given by the measured data on S. At time t = the solution of this problem 
will be equal to f{x). 

In even dimensions or when the sound speed is variable, the Huygens' principle 
does not hold anymore. Nevertheless, we could still try to rewind the wave in hopes of 
approximating the initial pressure. This is the main idea of the time reversal method 
and we will now give a precise description of the method in the general case of a 
variable sound speed and dimension n. 

To reconstruct the initial pressure f{x) inside B, we try to reverse the time (starting 
from time t = T and going back to t = 0) and solve the following problem: 

utt = c^(x)Au in B X [0,T] 

^ u{x,T)=p{x,T) 

Ut{x,T)=pt{x,T) 

. u\s{x,t) = g{x,t) on5x [0,T], 

where g = p\sx[o,T] is the restriction of p onto S x [0, T]. At time t = the solution 
of this problem equals f{xY. Obviously, this reconstruction requires the knowledge of 

'Note that u{x, t) is simply the restriction of p{x, t) onto B x [0, T]. 
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the solution to (1) on the cylinder 5* x [0, T] and inside B at time T. The values of p 
on the cylinder can be obtained from the measurements of the transducers placed on S. 
However, the values of the pressure inside B at time T are not known. Nevertheless, 
due to Theorem 2, we do know that the solution inside B decays with time. Therefore, 
after some time T > Tq it is reasonable to approximate p{x, T) and pt {x, T) with zero. 
This is how we arrive to the approximate reconstruction problem: 



Here {t) is a smooth cut-off function that equals to 1 in (— oo, T — e] (where T — e> 
To) and vanishes for t >T. As it will be explained below, at time f = the solution to 
(6) approximates f{x). In Section 3 we estimate the decay of this error with respect to 
the cut-off time T in the case of a non-trapping sound speed. 

3 An error estimate 

Consider the error e{x, t; T) = u{x, t) — v{x, t) that results from replacing the true 
solution of (5), u{x, t), with the time reversal solution (a; , t). Recall that w (a; , t) solves 
equation (6). Then, the error satisfies the following equation: 



To make the role of e more clear we will specify our choice of a cut-off function 
feit). Let ip{t) be a smooth function that equals to 1 on (— cxo, —1] and vanishes on 
[0, oo). When e < 1, we set (pe{t) = fii't- ~ T)/e) (Figure 1), and when e > 1, we 
choose Ve(i) = = This can also be written as (pe{t) = (p{{t — T)/a), 

where a = mm{e, 1}. 



' vtt = c^{x)Av kiBx [0,T] 

v{x,T) = 
Vt{x,T) = 

v\s{x,t) = g{x,t)^e{t) onSx[0,T], 



(6) 



< 



ett{x, t; T) = c'^{x)Ae{x, t; T) in B x [0, T] 

e{x,T;T) = p{x,T) 
et{x,T;T)=pt{x,T) 
^ e\s{x,t;T)=g,{x,t): = {1 - ipe{t))g{x,t) 



in 5 X [0,T]. 



(7) 
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Figure 1: A sketch of a cut-off function ^e{t)- 
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Let the sound speed c{x) satisfy the conditions described in Section 2.1, namely, 
c{x) > c > 0, c{x) G C°°(IR"), c{x) — 1 has a compact support and c(x) is non- 
trapping. Suppose also that the initial pressure f{x) belongs to L^(M") and is com- 
pactly supported (in B). Then the following theorem, which is the main result of this 
article, holds. 

Theorem 3. There exists Tq such that for any T > Tq and e > satisfying T — e > Tq, 
the error e{x, t; T) can be estimated as follows: 

• for even dimensions n 



Here C{e) = Cj min{e^, l},for some constant C depending on B and c{x). 
In particular, 



|e(0;T)||Hi(s) = \\f - vmmiB) <C{e){T - e)-^+^\\f\\L.,forevenn; 
||e(0;T)||^i(B) = ||/-^;(0)||hi(b) < C{e)e-'^'^-^^\\f\\L.(B),foroddn, 



where v{Q, x) is the approximation to f{x) obtained by time reversal with cut-off time 
T. 

Proof. Let E: J7^(S') , s > 0, be the operator of harmonic 

extension that produces a harmonic function E(j) in B from the Dirichlet boundary data 
(j) (e.g. [25]). Then w: = e — Eg^ satisfies 



As noted in Section 2.1, there exists a time Tq after which the solution p{x, t) to (1) is 
infinitely smooth in B x [T, oo) for any T > Tq. Let us choose such T and let e be 
such that T — e > Tq. Then the right hand side and the initial conditions of (8) are 
infinitely smooth functions. Indeed, p{x, T) is smooth, as we discussed above. Recall 
that is a smooth cut-off of the trace ofp{x, t). Moreover, it is supported in [T — e, T] 
- a time interval where p{x, t) has already become infinitely smooth. Thus Eg^ix, t) 
is a smooth function as well. 

Next, we will apply the following theorem to the solution of (8). 



Theorem 4. (e.g.[26]) Assume that G,H € C°°{B), F e C°°{B x [0, T]) and con- 
sider the initial boundary value problem 




wtt - (x) A«; = -Sffett in B X [0, T] 

w{x, T) = p{x, T) - Eg,{x, T) = p{x, T) - Eg{x, T) 
Wt{x,T) =pt{x,T) - Egst{x,T) =pt{x,T) - Egt{x,T) 
w\s{x,t) = onS'x[0,T], 



(8) 




Utt-c^{x)AU = F 
U = G, Ut = H 



in B X (0, T\ 
onB x{t = 0}, 
onS X [0, T] 



(9) 
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Here c{x) G C°° {B) is such that c{x) > 9 for some 6 > 0. Let the following compati- 
bility conditions hold for / = 1, 2, ... . 

Go := G e Hl{B), Hi :^ H e H^{B), 

G21 := ^^^(-,0) +c2(x)AG2;_2 G H^{B), (10) 
H21+1 ■■= ^^{■,0)+c^{x)AH2i-i e H^iB). 

Then (9) has a unique solution U G C°°{B x [0, T]) and 

^^^^i\WmH-(B) + ||f/t(i)IU^(B)) < 

Ce^'^{\\F\\L^O,T;L-{B)) + \\G\\hI(B) + \\H\\l-(b)), 
where the constants C and C\ depend on B and c{x). 

We will verify that the conditions of Theorem 4 are satisfied by (8), noting that the 

compatibility conditions must hold ait = T, since the direction of time is reversed. 
Indeed, as we already discussed, F = —Eg^^f, G = p(x, T) — Eg{x, T) and H = 
Pt{x,T) — Egt{x,T) are infinitely smooth. It is clear that Gq = G and Hi = H 
belong to H^^{B), as g{x,t) has been defined to be the trace of p{x,t). The higher 
order compatibility conditions hold as well: 



G2i{x) =-i^Eg,t^ix, T) + c2 AG2,-2 = -Ej^gix, T) + ^^€2(1-1) 
E^g{x,T) + (c2A)'Go = -E^g{x,T) + {c'' AYp{x,T). 



We used here that Ai? = 0. 

As we already know that p{x, T) is smooth in a neighborhood of T, we conclude 
that 

Mx,T) + {c^Ayp{x,T) 



G2iix)\s = 



= 0, 

s 



so G2; € i?o (-B) for / = 1, 2... Similarly, one can check that H21+1 € H^{B) for I = 
1,2.... 

Before applying Theorem 4 to the solution of (8) we notice that the energy in B, de- 
fined by ff(t) := i /g ^|Vw|^ + rr^{x) \wt\'^^ dx, stays constantin [0,T — a], where 
a = min{£, 1}. Indeed, using (8) one easily shows that 

S'{t) = / (Vw ■'Vwt + c~'^WtWtt) dx = - wtC~'^Egstt dx. 

J B J B 

The last integral vanishes in [0 , T — a] , because (?e = in [0 , T — a] and therefore, 
Eg^^^ = in this interval. Then, by Theorem (4), it follows that 



c{\\Eg + ||p(.,r)- Egi;T)f^,^^^ 

+ \\pt{;T)-Egt{;T)\\l,^j,^}, 
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where C depends on c{x) and B only, but not on T or e. Here and in what follows C 
will denote various constants, all of them independent on T or e. 
Asw = e — Eqs, the above inequalities imply 



pmax^dleWllffMB) + \\etit)\\mB)) < 

C m^^^l \\EgS)\\H-iB) + \\Eg,M\L^B)^ < 



^ "^^f^JW^aS^HHB) + \\Eg^t{t)\\mB)) + \\Eg,^^\\L2^T-a,T;LHB)) 

+ \\p{;T) - Eg{.,T)\\H^^^s) + \\Pt{;T) - i^^tl', 

(11) 

We used here that Eg^ = in [0, T - a]. 

LetT: H^{B) iJ*~^/^(S'), s > 1/2, denote the trace operator, which is known 
to be continuous (e.g.[26]). Then, using the continuity of the linear operators E and T, 
we can estimate the terms in (1 1) as follows. 

\\EgS)\\miB) < C\\ge{t)\\HU^s) = C\\{1 - Mt)mt)\\m/^s) < 

cm - MtmmHHB) < cMmH^By 

Similarly, 

\\Eg,M\L-iB) < \\Eg,,{t)\\Hi(B) < Cm - ^eM*)]* llffMs) < 

^{\\pmHHB) + \\PtmHHB)}, 

\\Eg,,,{t)\\mB) < \\Eg,,M\miB)<Cm-ip,)p]u\\miB) < 

^{MmHHB) + wPtmHHB) + wPttimHHB)} , 

\\Eg{T)\\HHB) < C\\p{T)\\HiiB) and \\Egt{T)h2iB) < C\\pt{T)\\HHB). 
Here we used that max [(^^(i)! < C/a and max < C" jo? for some constants 

C and C" . Then, the error can be estimated by 

„max^{||e(f)||j^i(B) + ||et(f)||i2(B)} < 

_max^y (Ib(i)llHMB) + IbtWIInMs)) + 

l|p|U2(T-a,T;Hi(B)) + lift |U2(T-a,T;Jfi(S)) + \Ptt\L'^(T -<x,T;m(B)) + 



c 



(12) 

Theorem 2 allows us to estimate the quantity inside the braces in right-hand side of 
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above inequality by a factor of 

^^^^rr.i'^oit) + m{t)] + \\Vo\\LHT-a,T) + \\Vl\\L^T-a,T) 
-a<t<T 

+ \\m\\L-{T-a,T) + vo{T) + m{T)y\f\\mB) 

The functions 77fe (t) are monotonically decreasing, and 772 (i) < ?7i(i) < Vo{t) provided 
t>l. Therefore, the right hand side of (12) is less than 

C/a^ {r?o(T - a) + \MlHt-c.,t)} WIUhb)- 
Taking into account that 

we arrive to 

m^^^(||e(t)|l^i(B) + WetimLHB)) < C(e)7?o(T - e)\\f\\L^^B), 

where r]a{T - e) = (T - e)-"+i for even n, r]o{T - s) = e'^^'^'^^ for odd n, and 
C{e) = C/min{£2,l}. 
This proves the theorem. 

4 Numerics of errors 

In this section we compare the errors of the numerical reconstructions of several phan- 
toms to their estimates given in Theorem 3. The numerical simulations were done in 
2D media with variable sound speed. Both cases of non-trapping and trapping sound 
speed are considered. 

For the computations we used the rectilinear finite difference scheme (as in [12, 
16]), implemented in Matlab. For both simulation of the phantom data and reconstruc- 
tion, we approximated the boundary of the unit circle, S, by the set of grid points 
closest to S and lying within B. For the simulation of phantom data the following 
problem was solved 

Ptt - c2 Ap in I? X (0, To] 
p(x,0) = /(x) Pt{x,Q)=Q 

P\dD = 0, 

where D = [—a, a]^ was a square containing S and large enough to ensure that no 
reflections off its boundary would reach S for time Tq. The values of the solution on 

5 were recorded for all time steps. The time Tq, the domain D and the space and time 
step-sizes were adjusted depending on the situation. For the reconstruction part, the 
same equation and difference scheme (but a different mesh) were used, this time on 
the square [—1.2, 1.2]^ instead of on D. To obtain the values of the boundary data at 
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Figure 2: A phantom (top) and its reconstructions using various radial sound speeds. 
Profiles of radial sound speeds are shown in the left column. In the right column are 
the reconstructions of the phantom, obtained by using the corresponding sound speeds 
from the left colunm. The first sound speed (second row) is non-trapping. The second 
one is a "trapping crater" sound speed. The third speed is a paraboloid, which causes 
severe trapping. The white dotted circles represent the observation surface S. 

the grid points used for reconstruction we used nearest point approximation. In order 
to observe the behavior of the error of reconstruction with respect to the cut-off time 
T, multiple reconstructions from different cut-off times were made and the error was 
graphed on a logarithmic scale. In both forward simulation and reconstruction, the 
spatial step-size h varied case by case and was of the order of 10~^. This corresponds 
to using several hundred detectors on the boundary of the unit circle. 

On Figure 2 examples of reconstructions using the time reversal method are shown. 
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Figure 3: Axial profile of a radial non-trapping sound speed (top), density plot of 
phantoms (left column) and the corresponding errors of reconstruction as functions 
of the cut-off time T (right column). The plots of the errors are in logarithmic scale, i.e. 
the horizontal axis represents InT and the vertical axis represents \n{H^ error). The 
dotted lines show the linear regression interpolation of the error. The decay of the error 
in the first example (second row) is of the order of T~^-^^, in the second one (third 
row) it is 



In the first case a variable non-trapping speed was used and the reconstruction is of a 
high quality. In the second example, the radial sound speed equals |a;| in the annulus 
A = {a;|0.5 < \x\ < 1} and is constant elsewhere. This "trapping crater" speed traps 
all rays corresponding to bicharacteristics that start at points (xq, Co) such that xq € A 
and a;o -L Co- The trapping leads to the blurring of the radial sides of the phantom, 
since those singularities do not reach the observation surface. A discussion of this 
phenomenon can be found in [12]. In the third example, shown in Figure 2, the sound 
speed equals |a;p -1-0.1 inside the unit circle. This speed exhibits a more severe trapping 
than the "trapping crater" speed, which leads to a much stronger blurring. In this case, 
for any xq there exists a cone of directions Co> such that the ray corresponding to the 
bicharacteristic starting at xq, Co is trapped [27]. 

In the next two figures we present various sound speeds and analyse the errors of 
reconstruction as functions of the cut-off time T. These figures show the particular 
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Figure 4: Density plot of a non-radial non-trapping sound speed (top), density plot of 
phantoms (left column) and the corresponding errors of reconstruction as functions 
of the cut-off time T (right column). The plots of the errors are in logarithmic scale. 
The dotted lines show the linear regression interpolation of the error. The decay of 
the error in the first example (second row) is of the order of T~^-^^, in the second one 
(third row) it is T-^ *^^. 



sound speed, phantoms and the corresponding errors, but not actual reconstructions, as 
we have already seen examples of such in Figure 2. The plot of each error shows only 
values of the error after some minimal time, specific to each sound speed, after which 
it makes sense to apply the time reversal procedure. Namely, this minimal time is taken 
to be the time a wave needs to cross the unit circle. After the error has decreased to a 
very small value it levels off, which is due to the discretization of the model. Thus, the 
error plots do not show these values. In order to easily observe the behaviour of the 
error, the plots are made in a logarithmic scale and the corresponding linear regression 
interpolations are graphed. 

Figures 3 and 4 show two examples of variable non-trapping sound speeds. In 
both cases two different phantoms are considered and the error of reconstruction 
as a function of the cut-off time T is plotted. In Theorem 3 we estimated that, in two 
dimensions, the error behaves as T^^. The examples from Figures 3 and 4 suggest 
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Figure 5: Axial profile of a radial "trapping crater" sound speed (top), density plot of 
phantoms (left column) and the corresponding errors of reconstruction as functions 
of the cut-off time T (right column). The plots of the errors are in logarithmic scale. 

even faster decay for those particular sound speeds. In the case shown in Figure 3 the 
error decays as T~^-^. The decay of the error in the second example (Figure 4) is 
even faster: the plots suggest that it is of the order of T~^-^. The next two examples 
deal with trapping sound speeds. In these cases the behaviour of the error depends also 
on the initial condition, as there is no uniform local energy decay unless additional 
smoothness is assumed (e.g. [28, 29]). The errors shown in Figures 5 and 6 decay 
as T increases. 

5 Final Remarks and Conclusion 

• In this text, we have concentrated on the errors resulting from the time reversal 
approximation. One is also interested in the stability of the method with respect 
to errors in the data. This question is answered by the standard results on stabihty 
of the mixed problem for the wave equation. If n{x,t) is the error in the mea- 
sured data g{x, t), then a stabihty estimate should bound a norm of e{x, 0; T), 
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Figure 6: Axial profile of a radial parabolic trapping sound speed (top), density plot of 
a phantom (left column) and the corresponding error of reconstruction as a function 
of the cut-off time T (right column). The plot of the error is in logarithmic scale. 

where e{x, t; T) solves the following mixed value problem: 



Since c{x) is assumed to be smooth, classical results (e.g. [30, Ch. 5, Sec. 5]) 
give, for instance, the following estimate: 



l|e(a;,0;T) 11^2(3) < C||n(x,t)||^i+^.i+^([o,T]xS) 
Here H^+^^^+^ {[0,T] x S) ^ L^{0,T; H^+^{S)) n H^+^iO,T: 1^(3)) and 



6 > 0. In a nutshell, the reconstruction is as stable as Radon inversion in the 
corresponding dimension. 

• Our numerical examples show decay of the error also in the trapping speed 

case. One could probably get some estimates on this decay, assuming sufficient 
smoothness of the function f{x) to be reconstructed and using the results of [29]. 
It is unlikely that one can get such estimates without a smoothness assumption. 

To conclude, in this paper we have provided error estimates of the time reversal justi- 
fying it in any dimension under non-trapping condition. The numerical examples agree 
with the error estimates in the cases of non-trapping speeds and show that the time 
reversal works even when the speed is trapping. 



' ett{x, t; T) = c'^{x)Ae{x, t;T) in B x [0, T] 

e{x,T;T) = 

et{x,T;T)=0 
. e\s{x,t;T) =ne{x,t): = {1 - (ps{t))n{x,t) 



inSx [0,T]. 



(13) 
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